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We propose a new Monte Carlo algorithm for the numerical study of general lattice models in 
Hamiltonian form. The algorithm is based on an initial Ansatz for the ground state wave function 
depending on a set of free parameters which are dynamically optimized. The method is discussed 
in details and results are reported from explicit simulations of U(l) lattice gauge theory in 1 + 1 
dimensions. 
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I. INTRODUCTION 

The lattice formulation of Quantum Field Theory is a 
rigorous theoretical framework for the study of non per- 
turbativc phenomena and its applications to phenomcno- 
logically relevant models, like QCD, are expected to play 
a major role in the next future. In most numerical 
studies, the lattice version of a given model is built in 
the Lagrangian formulation pi. However, for certain 
"static" problems, like hadron spectroscopy, the Hamil- 
tonian approach is more suited because time is kept 
as continuous as possible and the spectrum is deformed 
in a minimal way. Actually, even the most accurate La- 
grangian calculations ||] evidence the advantages of using 
anisotropic lattices with smaller spacings in the temporal 
direction, a typical feature of the Hamiltonian approach. 

Another advantage of the Hamiltonian formulation is 
the possibility of easily exploiting any knowledge about 
the ground state wave function \I/o i n order to improve 
the performance of numerical calculations. In partic- 
ular, in this Letter, we are concerned with the popu- 
lar and all purpose Green Function Monte Carlo algo- 
rithm (GFMC) @ where an approximate version of ^q, 
the trial wave function ^>t, allows to speed up the al- 
gorithm convergence by means of Importance Sampling 
techniques K). 

In a realistic problem, ^Bx depends on several free pa- 
rameters a = (ai, . . . , a at) that parametrize interaction 
terms responsible for correlations expected to show up in 
the ground state. The development of methods for the 
optimization of the free parameters a is a key problem 
and an active field of research Q . The simplest approach 
is to perform a variational calculation of the mean energy 
E v (a) = (* T (a)|/f|*T(a)) where H is the Hamilto- 
nian. The quantity Ey (a) is then minimized with respect 
to a. Since each evaluation of Ey (a) requires a separate 
Monte Carlo calculation, this is an expensive procedure 
for which special tricks have been devised (e.g. correlated 
sampling ) . A particularly interesting approach is that 
of M where the authors propose an algorithm for the au- 
tomatic dynamical optimization of the free variational 



parameters a. 

In this Letter, we abandon the variational approxi- 
mation and consider instead a full Monte Carlo calcu- 
lation for which we propose a new strategy to optimize 
a. We describe an adaptive algorithm which converges 
on-line to an optimal set of parameters a* which min- 
imizes the statistical error of the full ( not variational) 
Monte Carlo simulation. The non-adaptive core of the 
algorithm is a GFMC algorithm with Importance Sam- 
pling and Stochastic Reconfiguration ||. 

Let us consider a quantum mechanical point particle 
with Hamiltonian 



H = H + V(q), 
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Ho = ^P 2 



(1) 



where q = (q±, . . . ,qd) is the position in R d and p — 
(pi, . . . ,Pd) is the associated momentum satisfying the 
canonical commutation rules [qi,qj] = = [pi,pj] and 
[9»>2>j] = iSij- Let *T(q, a) = expF(q, a) be a positive 
approximation of the ground state wave function depend- 
ing on some parameters a = [ai, . . . , ajy). We perform 
a unitary transformation on H based on and build 
the new isospectral Hamiltonian 



H = V T H 
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(2) 



The Feynman propagator D = (q"| exp(— £_ff)|q') associ- 
ated with H can be written 



D(q",q',a, £ ) = 



1 



1 



(27re) rf / 2 Vl + eV 2 ^ 
x exp {-(xifciSij + edfjF^Xj - ev} + 0(e 2 ), 

Xi = Qi zdiF, 
where 



(3) 
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1. 



V(q, a) = V(q) - -(V q F(q, a)) 2 - -V 2 F(q, a). (4) 

By standard arguments, we conclude that D may be ex- 
pressed in terms of a weighted average 
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5(q",q',a,t) = / Vq(t) e £ ^'), 

Jq(0)=q',q(t)=q" 

(5) 

where I? is the formal limiting measure (when e — > 0) 
defined on the stochastic paths generated according to 
the Langevin equation 



ediF(q n ,a) 
1 
2 



(6) 



e5?-F(q n! a)z„ ;i ), 



where <\ n i denotes the i-th component of q„ and z ra 
are uncorrelated Gaussian random points in R d with 
unit variance. The GFMC algorithm is expressed in 
a concise and formal way by Eq. , but the ac- 
tual calculation of the weighted averages must keep un- 
der control the variance of the path-dependent weights 
exp{— J Q drV(q(r), a)} which explodes as t — > +oo. An 
efficient technique to solve this problem is Stochastic Re- 
configuration |j] which implements a kill and branch se- 
lection on the paths with the desirable feature of dealing 
always with a fixed size population of walkers. To this 
aim, a finite collection of K walkers, an ensemble, is in- 
troduced 

£ = {(q (n) (t),o;W(t))} 1 <„< if , (7) 
where the weights uA n ' are defined by 

uj {n \t) =exp- / drF(qW(r),a). (8) 



The variance of the weights over the ensemble £ is 



W(t) = Var uj(t) 



1 K 
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fc=l 



K "52 



,(fc) 
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fc=l 



and the average of a function /(q) over £ is 



(9) 



(10) 



fc=l 



When W(i) becomes too large, £ is transformed into a 
new ensemble £' with zero variance W and the same av- 
erages, at least in the K — > oo limit. Simulations can be 
extended to arbitrarily large times with the drawback of 
a systematic error vanishing with K and an extrapolation 
to K — > oo is required. More details on this procedure 
can be found in 0. 

Statistical errors in such a Monte Carlo simulation are 
related to the fluctuations of V r (q(r),a). To see this, 
let us consider for instance the calculation of the ground 
state energy Eq. A simple algorithm estimates the limit 



Eq = lim t ^ +oc E(t, q', a) where, for arbitrary q' and a , 
we define 



£(t,q',a) = --log / D(q",q',a,*)dq" = (F(q(t),a)>, 

(11) 

((■) is the average over the weighted Langevin trajecto- 
ries as in Eq. (||)). If is an exact eigenstate of H 
with eigenvalue E, then we obtain V = E and the above 
formula gives E(t) — E with zero variance, namely no 
statistical error. 

In the more general case of non optimal "fy 
a simulation performed with a sufficiently small e and a 
population of K walkers will provide after S Monte Carlo 
steps only an approximate estimator Eq(S, K, a) of Eq, 
that is a random variable with the asymptotic properties 

(E Q (S,K, a )) =E Q + ^ + o(K- a ), a>0, (12) 



Var E (S, K,a) 



C2Cg>) 

Vs ' 



where averages are over Monte Carlo realizations. The 
average of Eq extrapolated at K — > oo is exact and inde- 
pendent on the trial parameters a whereas C2 (K, a) is in 
general strongly dependent on them and is expected to 
vanish when K — * oo. 

If the family ^t((1, a) includes the exact ground state 
at the special point a = a*, then we know that Ci(a*) = 
C2(K, a*) = 0; in a less optimal situation, motivated 
by [fL0|, we seek a minimum of C2. 

To establish an adaptive algorithm, we let a — > {a„} 
be a dynamical parameter of the simulation and propose 
to update it together with the ensemble according to the 
equations 

q«+M = Qn,i + ediF(q n , a„) + 



+ y/e(x n ,i + -£9|-F(q n ,a„)z n j), 
a„+i = a„ - ?7V a JF n (a n ), 
T n (si n ) = Var £ V"(q„,a„), 



(13) 



where 7/ is a constant parameter and a n is the value of 
a at the n-th update. In other words, we implement a 
local minimization of the weight variance as a driving 
mechanism for the free parameters. 

The coupled set of equations Eqs. ([l3|) for the evolution 
of a and the random walkers is non linear and discrete. 
The optimization of a and the GFMC are thus linked 
together. This procedure can be successfully checked in 
trivial quantum mechanical examples. Here we discuss a 
non trivial application to test convergence and stability. 
The model we study is U(l) lattice gauge theory in 1 + 1 
dimensions. Following the notation of [jllj (see also C^] 
for other applications of GFMC to lattice gauge theory) , 
the Hamiltonian of the model is 
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(14) 



where (3 is the coupling constant, L is the spatial lat- 
tice size, are link phases around the p-th plaque- 
tte and 4>p — @i,p + @2.p+i — #3, p — ®2,p is the plaquette 
gauge invariant angle. An accurate variational estimate 
of the ground state energy per plaquette Eq/L is ob- 
tained with the gauge invariant Ansatz Wo(^i, . . . , 4>l) — 
exp ^A J2p=i cos ^p) leading to 



E, 



(var) 



L 



= 0- 
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(A-/3 2 ) 



h(2X) 
/o(2A) 



(15) 



where /„ is the n-th modified Bessel function. 

We choose a rather general gauge invariant trial wave 
function — exp F with inter plaquette correlations of 
the following form 



F = cos cj) p + a-2 cos20 p + (16) 

2 

+ £(°3,fc C0S (^P + ^P+k) + 04, k COS(0 P - <f>p+k))}, 
k=l 

which depends on 6 free parameters a,\, &2, &3,i, 0-3,2, 
04,1, 04,2- The particular case 03,2 = 04,2 = with 
only on-site and nearest-neighbour terms is discussed 
in |flj where the 0{(3 A ) perturbative optimal values are 

pert 

given: a\ 



/? 2 /2, af rt = -/? 4 /32, afff = /3 4 /24 and 
a pert _ fog g roun( j state energy is measured 

by using only ai which is optimized by trials and errors. 
In this Letter, we shall discuss the non perturbative be- 
haviour of the full six parameters at the equilibrium point 
reached automatically by the algorithm. 

We perform simulations with ensembles of K = 10, 30 
and 50 walkers on a system with spatial lattice size L = 8 
and time step e = 0.015. In the Lagrangian formulation 
this would correspond to a simulation on a lattice with a 
very large extension in the temporal direction. The cou- 
pling (3 is varied through the values 0.5, 1.0, 1.5, 2.0 and 
2.5. We execute stochastic reconfiguration of the ensem- 
ble each r(/3) temporal steps with r(/3) ranging between 
a maximum 40 at [3 — 0.5 and a minimum 20 at (3 = 2.5. 
The learning parameter is r\ — 0.0005 and all the free 
parameters are zero at the first iteration. The number 
of Monte Carlo iterations is around 5 x 10 4 depending 
slightly on (3. 

In Fig. (1) we show the Monte Carlo time history of the 
free parameters a for the run with (3 = 0.5 and K = 30. 
After about 300 iterations, the parameters reach an op- 
timum value around which they oscillate with small fluc- 
tuations. If these fluctuations were too noisy, they could 
be removed by Stochastic Gradient Approximation tech- 
niques, namely by letting r\ be a time dependent posi- 
tive sequence {rj n } n >o vanishing with n — > 00 under the 



constraints J^ n V = 00 an d V 2 < 00 @] ■ The corre- 
sponding energy time history is shown in Fig. (2) where 
the fluctuations of the local measurements are shown to 
decrease very rapidly. 

The optimal values of the six parameters as functions 
of (3 at the largest K used are shown in Tab. (||). The 
leading perturbative expansion of a can be checked to 
be correct only at the smallest (3 and overestimates a* at 
larger (3. The next-to-nearest neighbour terms are rather 
small being at (3 — 2.5 roughly 1% of the nearest neigh- 
bour ones. The qualitative picture discussed in Q is 
confirmed with dominating plaquette anticorrelation in 
the ground state and next-to-nearest neighbour effects 
below the percent level. Here, we stress again, the de- 
termination of the optimal set of parameters a* is com- 
pletely automatic. 

About the ground state energy, we show in Tab. (0) 
the estimates computed with three values of K to show 
the very small systematic residual error associated with 
the finite size walker population. We also show the Monte 
Carlo results from jll| with which we agree within errors 
as well as the variational bound. It is remarkable that 
a good estimate of the ground state energy is obtained 
with the very small number of walkers K = 10. This can 
be intepreted as a signal that for this admittedly sim- 
ple model the proposed six parameters wave function is 
rather accurate. 

To conclude, let us remark that the above method is 
rather general and is applicable to dynamically optimize 
the free parameters of a many body trial wave func- 
tion for any model that can be studied by GFMC with 
Stochastic Reconfiguration. This includes for instance 
pure gauge SU(2) and SU(3) lattice gauge theory in any 
dimension. 
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FIG. 1. Monte Carlo time history of the free parameters a 
for the run with /3 = 0.5 and K = 30. 
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FIG. 2. Monte Carlo time history of the energy measure- 
ments during the first steps of the run with (5 — 0.5 and 
K = 30. 
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TABLE I. Optimal value of the 6 free parameters 01,02,03,1,03,2,04,1,04,2 and corresponding minimum standard deviation 
cr_B of the ground state energy estimator. We show the numerical results obtained with 50 walkers. 



p 


ai 


02 


03,1 


a.3,2 


04,1 


04,2 


OB 


0.5 


0.1245271(7) 


-0.0019590(2) 


0.0024544(4) 


-0.0000374(3) 


-0.0015535(3) 


-0.0000609(3) 


0.0002 


1.0 


0.46287(4) 


-0.022487(6) 


0.02597(1) 


0.00074(1) 


-0.015262(8) 


-0.00108(1) 


0.003 


1.5 


0.8327(2) 


-0.05406(3) 


0.06255(5) 


0.00319(4) 


-0.02999(4) 


-0.00414(4) 


0.009 


2.0 


1.1541(6) 


-0.07794(5) 


0.0997(1) 


0.00618(7) 


-0.03921(7) 


-0.00866(6) 


0.008 


2.5 


1.4447(5) 


-0.09663(7) 


0.1370(1) 


0.01065(9) 


-0.0459(1) 


-0.0128(1) 


0.008 



TABLE II. Comparison among the different estimates (the proposed method, standard Monte Carlo and variational) of the 
ground state energy. 





77iAdaptive 
- £i MC.K = 10 


^Adaptive 
-^^10^ = 30 


^Adaptive 
-^MC.K^SO 


-£<MC, Barnes et al. 


-^Variational 


0.5 


0.468930(2) 


0.468932(1) 


0.4689327(8) 


0.4690(1) 


0.4690 


1.0 


0.77032(3) 


0.77020(2) 


0.77024(2) 


0.7697(2) 


0.7746 


1.5 


0.88226(8) 


0.88209(5) 


0.88188(4) 


0.8823(8) 


0.9005 


2.0 


0.91718(9) 


0.91694(5) 


0.91705(4) 


0.916(1) 


0.9435 


2.5 


0.93245(8) 


0.93242(5) 


0.93230(3) 


0.929(2) 


0.9594 
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